subroutine CalcJacobians()
use DataJacobians
use Data_cfoam
use Data_Mesh, only: nFaces, nInternalFaces, meshFaceOwner,meshFaceNeighbour, &
                     opposingOwnerFaceLabel,opposingNeighbourFaceLabel, &
                     meshCells,meshFaces,meshPoints
implicit none

    integer :: iFace, iFaceB,iFaceF, iFlocal
    integer :: iFaceBbeta(2), iFaceFbeta(2)
    integer :: iFaceBgamma(2),iFaceFgamma(2)
    
    integer :: iCellB,iCellF
    integer :: i,j,k, n, a,b, iPoint1,iPoint2
    
    real(8) :: Lalfa(3,2),Lbeta(3,2),Lgamma(3,2)
    real(8) :: det
    
    logical :: OK_FLAG

    allocate(Jtb(3,3,nFaces))
    allocate(Jtcb(3,3,nFaces))
    allocate(Jtn(3,3,nFaces))
    allocate(Jtcf(3,3,nFaces))
    allocate(Jtf(3,3,nFaces))

    allocate(invJtb(3,3,nFaces))
    allocate(invJtcb(3,3,nFaces))
    allocate(invJtn(3,3,nFaces))
    allocate(invJtcf(3,3,nFaces))
    allocate(invJtf(3,3,nFaces))

    allocate(invJtbInvJb(3,3,nFaces))
    allocate(invJtcbInvJcb(3,3,nFaces))
    allocate(invJtnInvJn(3,3,nFaces))
    allocate(invJtcfInvJcf(3,3,nFaces))
    allocate(invJtfInvJf(3,3,nFaces))

    do iFace = 1,nInternalFaces
        iCellB = meshFaceOwner(iFace)
        iCellF = meshFaceNeighbour(iFace)

        iFaceB = opposingOwnerFaceLabel(iFace)
        iFaceF = opposingNeighbourFaceLabel(iFace)
        ! alfa-direction = iCellB - iFace - iCellF


! ********************************************************************************* !

        ! Def beta and gamma directions in cell B
        iFaceBbeta(1:2) = 0
        do iFlocal = 1,6
            i = meshCells(iFlocal,iCellB) 
            if(i.ne.iFace.and.i.ne.iFaceB)then
                iFaceBbeta(1) = i
                exit
            end if
        end do
        do iFlocal = 1,6
            i = meshCells(iFlocal,iCellB) 
            if(i.ne.iFace.and.i.ne.iFaceB.and.i.ne.iFaceBbeta(1))then
                n = 0
                do a = 1,4
                    iPoint1 = meshFaces(a,i)
                    do b = 1,4
                        iPoint2 = meshFaces(b,iFaceBbeta(1))
                        if(iPoint1.eq.iPoint2)then
                            n = n + 1
                        end if
                    end do
                end do

                if(n.eq.0)then
                    iFaceBbeta(2) = i
                    exit
                end if
            end if
        end do
        
        if(iFaceBbeta(1).eq.0.or.iFaceBbeta(2).eq.0)then
            Stop '@@FATAL ERROR: iFaceBbeta has not been defined!'
        end if
        
        iFaceBgamma(1:2) = 0
        k=1
        do iFlocal = 1,6
            i = meshCells(iFlocal,iCellB) 
            if(i.ne.iFace.and.i.ne.iFaceB.and.i.ne.iFaceBbeta(1).and.i.ne.iFaceBbeta(2))then
                iFaceBgamma(k) = i
                if(k.eq.2) exit
                k = k + 1
            end if
        end do

        if(iFaceBgamma(1).eq.0.or.iFaceBgamma(2).eq.0)then
            Stop '@@FATAL ERROR: iFaceBgamma has not been defined!'
        end if

        ! Def beta and gamma directions in cell F
        do i = 1,2
            call CalcJoinedFace(iCellF,iFace,iFaceBbeta(i),iFaceFbeta(i))
            call CalcJoinedFace(iCellF,iFace,iFaceBgamma(i),iFaceFgamma(i))
        end do

! ********************************************************************************* !

        ! Def Jacobians
        ! for Bc
        call FaceCenterDef(iFace,Lalfa(1:3,2))
        call FaceCenterDef(iFaceB,Lalfa(1:3,1))
        
        do i = 1,2
            call FaceCenterDef(iFaceBbeta(i),Lbeta(1:3,i))
            call FaceCenterDef(iFaceBgamma(i),Lgamma(1:3,i))
        end do
                
        call CalcJacobian(Lalfa,Lbeta,Lgamma,Jtcb(1:3,1:3,iFace))
        
        ! for Fc
        call FaceCenterDef(iFaceF,Lalfa(1:3,2))
        call FaceCenterDef(iFace,Lalfa(1:3,1))
        
        do i = 1,2
            call FaceCenterDef(iFaceFbeta(i),Lbeta(1:3,i))
            call FaceCenterDef(iFaceFgamma(i),Lgamma(1:3,i))
        end do
                
        call CalcJacobian(Lalfa,Lbeta,Lgamma,Jtcf(1:3,1:3,iFace))

        ! for b
        !call CellCenterDef(iCellB,Lalfa(1:3,2))
        call FaceCenterDef(iFace,Lalfa(1:3,2))
        call FaceCenterDef(iFaceB,Lalfa(1:3,1))
        
        do i = 1,2
            call CalcEdgeCenter(iFaceB,iFaceBbeta(i),Lbeta(1:3,i))
            call CalcEdgeCenter(iFaceB,iFaceBgamma(i),Lgamma(1:3,i))
        end do

        call CalcJacobian(Lalfa,Lbeta,Lgamma,Jtb(1:3,1:3,iFace))
        
        ! for f
        !call CellCenterDef(iCellF,Lalfa(1:3,2))
        call FaceCenterDef(iFace,Lalfa(1:3,1))
        call FaceCenterDef(iFaceF,Lalfa(1:3,2))
        do i = 1,2
            call CalcEdgeCenter(iFaceF,iFaceFbeta(i),Lbeta(1:3,i))
            call CalcEdgeCenter(iFaceF,iFaceFgamma(i),Lgamma(1:3,i))
        end do

        call CalcJacobian(Lalfa,Lbeta,Lgamma,Jtf(1:3,1:3,iFace))
        
        ! for n
        call CellCenterDef(iCellF,Lalfa(1:3,2))
        call CellCenterDef(iCellB,Lalfa(1:3,1))
        do i = 1,2
            call CalcEdgeCenter(iFace,iFaceBbeta(i),Lbeta(1:3,i))
            call CalcEdgeCenter(iFace,iFaceBgamma(i),Lgamma(1:3,i))
        end do

        call CalcJacobian(Lalfa,Lbeta,Lgamma,Jtn(1:3,1:3,iFace))
        
        !call CalcInvMatrix33(Jtn(1:3,1:3,iFace),invJtn(1:3,1:3,iFace))
        
! ********************************************************************************* !
        ! Definition inverse Jacobians

        ! invJtb
        call M33INV(Jtb(1:3,1:3,iFace),invJtb(1:3,1:3,iFace), OK_FLAG)
        if(OK_FLAG.eqv..false.)then
            stop '@@FATAL ERROR: Inverse matrix 3 by 3 has not been calculated'
        end if

        ! invJtcb
        call M33INV(Jtcb(1:3,1:3,iFace),invJtcb(1:3,1:3,iFace), OK_FLAG)
        if(OK_FLAG.eqv..false.)then
            stop '@@FATAL ERROR: Inverse matrix 3 by 3 has not been calculated'
        end if
        
        ! invJtn
        call M33INV(Jtn(1:3,1:3,iFace),invJtn(1:3,1:3,iFace), OK_FLAG)
        if(OK_FLAG.eqv..false.)then
            stop '@@FATAL ERROR: Inverse matrix 3 by 3 has not been calculated'
        end if

        ! invJtcf
        call M33INV(Jtcf(1:3,1:3,iFace),invJtcf(1:3,1:3,iFace), OK_FLAG)
        if(OK_FLAG.eqv..false.)then
            stop '@@FATAL ERROR: Inverse matrix 3 by 3 has not been calculated'
        end if

        ! invJtf
        call M33INV(Jtf(1:3,1:3,iFace),invJtf(1:3,1:3,iFace), OK_FLAG)
        if(OK_FLAG.eqv..false.)then
            stop '@@FATAL ERROR: Inverse matrix 3 by 3 has not been calculated'
        end if

! ********************************************************************************* !
        ! Defintion products JtInvJ
    
        ! invJtbJb
        call CalcProductAAt(invJtb(1:3,1:3,iFace),invJtbInvJb(1:3,1:3,iFace))
        
        ! invJtcbJcb
        call CalcProductAAt(invJtcb(1:3,1:3,iFace),invJtcbInvJcb(1:3,1:3,iFace))
        
        ! invJtnJn
        call CalcProductAAt(invJtn(1:3,1:3,iFace),invJtnInvJn(1:3,1:3,iFace))
        
        ! invJtcfJcf
        call CalcProductAAt(invJtcf(1:3,1:3,iFace),invJtcfInvJcf(1:3,1:3,iFace))
        
        ! invJtfJf
        call CalcProductAAt(invJtf(1:3,1:3,iFace),invJtfInvJf(1:3,1:3,iFace))
        
! ********************************************************************************* !        
    end do

end subroutine CalcJacobians

subroutine CalcProductAAt(A,AAt)
implicit none
    real(8),intent(in) :: A(3,3)
    real(8),intent(out) :: AAt(3,3)
    
    integer :: i,j,k
    
    AAt(1:3,1:3) = 0.0d0
    do j = 1,3
        do i = 1,3
            do k = 1,3
                AAt(i,j) = AAt(i,j) + A(i,k)*A(j,k)
            end do
        end do
    end do
    
end subroutine CalcProductAAt

! Inverse matrix 33 calculations
      SUBROUTINE M33INV (A, AINV, OK_FLAG)

      IMPLICIT NONE

      real(8), DIMENSION(3,3), INTENT(IN)  :: A
      real(8), DIMENSION(3,3), INTENT(OUT) :: AINV
      LOGICAL, INTENT(OUT) :: OK_FLAG

      real(8), PARAMETER :: EPS = 1.0D-10
      real(8) :: DET
      real(8), DIMENSION(3,3) :: COFACTOR


      DET =   A(1,1)*A(2,2)*A(3,3)  &
            - A(1,1)*A(2,3)*A(3,2)  &
            - A(1,2)*A(2,1)*A(3,3)  &
            + A(1,2)*A(2,3)*A(3,1)  &
            + A(1,3)*A(2,1)*A(3,2)  &
            - A(1,3)*A(2,2)*A(3,1)

      IF (ABS(DET) .LE. EPS) THEN
         AINV = 0.0D0
         OK_FLAG = .FALSE.
         RETURN
      END IF

      COFACTOR(1,1) = +(A(2,2)*A(3,3)-A(2,3)*A(3,2))
      COFACTOR(1,2) = -(A(2,1)*A(3,3)-A(2,3)*A(3,1))
      COFACTOR(1,3) = +(A(2,1)*A(3,2)-A(2,2)*A(3,1))
      COFACTOR(2,1) = -(A(1,2)*A(3,3)-A(1,3)*A(3,2))
      COFACTOR(2,2) = +(A(1,1)*A(3,3)-A(1,3)*A(3,1))
      COFACTOR(2,3) = -(A(1,1)*A(3,2)-A(1,2)*A(3,1))
      COFACTOR(3,1) = +(A(1,2)*A(2,3)-A(1,3)*A(2,2))
      COFACTOR(3,2) = -(A(1,1)*A(2,3)-A(1,3)*A(2,1))
      COFACTOR(3,3) = +(A(1,1)*A(2,2)-A(1,2)*A(2,1))

      AINV = TRANSPOSE(COFACTOR) / DET

      OK_FLAG = .TRUE.

      RETURN

      END SUBROUTINE M33INV


! Definition the center of the edge common for faces iFace1 and iFace2
subroutine CalcEdgeCenter(iFace1,iFace2,R)
use Data_Mesh, only: meshFaces, meshPoints
implicit none
    integer,intent(in) :: iFace1,iFace2
    real(8),intent(out) :: R(3)
    
    integer :: a,b, n, iPoints(2)
    integer :: iPoint1,iPoint2
    
    R(1:3) = 0.0d0

    n = 0
    do a = 1,4
        iPoint1 = meshFaces(a,iFace1)
        do b = 1,4
            iPoint2 = meshFaces(b,iFace2)
            if(iPoint1.eq.iPoint2)then
                n = n + 1
                if(n.gt.2)then
                    Stop '@FATAL ERROR: iPoints has not been found!'
                end if
                iPoints(n) = iPoint1
            end if
        end do
    end do

    do a = 1,2
        R(1:3) = meshPoints(1:3,iPoints(a))
    end do

end subroutine CalcEdgeCenter

subroutine CellCenterDef(iCell,R)
use Data_Mesh, only: meshPoints, meshCellPoints
implicit none
    integer,intent(in) :: iCell
    real(8),intent(out) :: R(3)
    
    integer :: i
    
    R(1:3) = 0.0d0
    do i = 1,8
        R(1:3) = R(1:3) + meshPoints(1:3,meshCellPoints(i,iCell))
    end do
    R(1:3) = 0.125d0 * R(1:3)

end subroutine CellCenterDef


subroutine FaceCenterDef(iFace,R)
use Data_Mesh, only: meshPoints, meshFaces
implicit none
    integer,intent(in) :: iFace
    real(8),intent(out) :: R(3)
    
    integer :: i
    
    R(1:3) = 0.0d0
    do i = 1,4
        R(1:3) = R(1:3) + meshPoints(1:3,meshFaces(i,iFace))
    end do
    R(1:3) = 0.25d0 * R(1:3)

end subroutine FaceCenterDef

! Jt = J !
!          / dx/da dy/da dz/da \  Lalfa
! J(3:3) = | dx/db dy/db dz/db |  Lbeta
!          \ dx/dc dy/dc dz/dc /  Lgamma
! J(2,1) = dy/da
! J(1,2) = dx/db
!
subroutine CalcJacobian(Lalfa,Lbeta,Lgamma,J)
implicit none
    real(8),intent(in)  :: Lalfa(3,2), Lbeta(3,2), Lgamma(3,2)
    real(8),intent(out) :: J(3,3)

    integer :: i
    
    do i = 1,3
        J(1,i) = Lalfa(i,2) - Lalfa(i,1)
        J(2,i) = Lbeta(i,2) - Lbeta(i,1)
        J(3,i) = Lgamma(i,2) - Lgamma(i,1)
    end do

end subroutine CalcJacobian

subroutine CalcJoinedFace(iCell,iFace,iFaceB,iFaceF)
use Data_Mesh, only: meshCells, meshFaces
implicit none
    integer,intent(in) :: iCell,iFace,iFaceB
    integer,intent(out) :: iFaceF

    integer :: i, a,b, n, jFace
    integer :: iPoint1, iPoint2

    iFaceF = 0

    do i = 1,6
        jFace = meshCells(i,iCell)
        if(jFace.eq.iFace)cycle
        
        n = 0
        do a = 1,4
            iPoint1 = meshFaces(a,jFace)
            do b = 1,4
                iPoint2 = meshFaces(b,iFaceB)
                if(iPoint1.eq.iPoint2)then
                    n = n + 1
                end if
            end do
        end do

        if(n.eq.2)then
            iFaceF = jFace
            return
        end if

    end do

    Stop '@FATAL ERROR: iFaceF has not been found!'

end subroutine CalcJoinedFace
